##########################################################################################

library(data.table)
library(dplyr)
library(ggplot2)
library(ggpubr)
library(optparse)
library(tidyr)

##########################################################################################

option_list <- list(
  make_option(c("--info_file"), type = "character"),
  make_option(c("--out_path"), type = "character")
)

###########################################################################################

parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

if(1!=1){
  info_file <- "/public/home/xxf2019/20220915_gastric_multiple/dna_combinePublic/public_ref/combine/MutationInfo.combine.addMolecularSubType.Race.tsv"
}

info_file <- opt$info_file
out_path <- opt$out_path

###########################################################################################

dat_info <- data.frame(fread(info_file))

###########################################################################################
## 构建All
dat_info_all <- dat_info
dat_info_all$From <- "All"

dat_info <- rbind(dat_info,dat_info_all)

##MSI和POLE合在一起
##unknown和EBV和在一起
dat_info$Molecular.subtype <- ifelse(dat_info$Molecular.subtype %in% c("MSI","POLE"),"MSI/POLE",
                                     ifelse(dat_info$Molecular.subtype %in% c("unknown","EBV"),"Unknown",dat_info$Molecular.subtype))

## IGC+DGC的混合型,分别拆分加在样本上面
#tmp <- subset(dat_info , Class == "IGC + DGC")
#tmp1 <- tmp
#tmp1$Class <- "IGC"
#tmp2 <- tmp
#tmp2$Class <- "DGC"
#dat_info <- rbind( dat_info , tmp1 , tmp2 )

###########################################################################################
## 按照患者计算
dat_plot <- dat_info %>%
  group_by(Class,From,Molecular.subtype) %>%
  summarise(count=n()) %>% 
  group_by(Class,From) %>%
  mutate(count_all=sum(count),
         ratio=count/count_all)

dat_plot <- subset(dat_plot,dat_plot$Class %in% c("IGC","DGC"))
dat_plot$p <- ""

for(i in unique(dat_plot$From)){
  dat_p <- subset(dat_plot,From==i)
  dat_p <- dat_p[,c("Class","Molecular.subtype","count")]
  matrix_data <- dat_p %>%
    group_by(Class) %>%
    spread(Molecular.subtype, count,fill = 0)
  if(nrow(matrix_data)==2){
    dat_plot$p[dat_plot$From==i][1] <- fisher.test( matrix_data[,-1] )$p.value
  } else {
    dat_plot$p[dat_plot$From==i] <- ""
  }
}

trans <- function(num){
  up <- floor(log10(num))
  down <- round(num*10^(-up),2)
  text <- paste("P == ",down," %*% 10","^",up)
  return(text)
}
dat_plot$p <- as.numeric(dat_plot$p)
dat_plot$value_text <- round(dat_plot$ratio , 2) * 100 
dat_plot$p_text <- ifelse( dat_plot$p < 0.01 , trans(dat_plot$p) , paste0( "P == " , round(as.numeric(dat_plot$p) , 3) )  )
dat_plot$p_text <- ifelse( is.na(dat_plot$p_text) , "" , dat_plot$p_text  )


###########################################################################################
col <- c(
  rgb(204,116,92,alpha=255,maxColorValue=255),
  rgb(62,62,93,alpha=255,maxColorValue=255),
  rgb(229,198,143,alpha=255,maxColorValue=255),
  "grey"
)

names(col) <- c("CIN","GS","MSI/POLE","Unknown")
dat_plot$Molecular.subtype <- factor(dat_plot$Molecular.subtype,levels=c("CIN","GS","MSI/POLE","Unknown"))
dat_plot$Class <- factor(dat_plot$Class,levels=c("IGC","DGC"))
dat_plot$Class_n <- paste0( dat_plot$Class , "\n(" , dat_plot$count_all , ")")
dat_plot$Class_n <- factor( dat_plot$Class_n , levels = unique(dat_plot$Class_n )[order(unique(dat_plot$Class_n ) , decreasing=T )] )

p <- ggplot(dat_plot , mapping = aes(Class_n , ratio , fill = Molecular.subtype)) +
  geom_bar(position = "stack", stat = "identity") + 
  theme_bw() +
  facet_grid(.~From,space='free_x',scales='free_x') +
  theme(strip.text.x = element_text(size = 7 , face = 'bold'))+
  theme(panel.grid=element_blank())+
  labs(y = 'Proportions (%)') +
  geom_text(aes(label=p_text , y = 1.05 ,x=1.5),size=3,parse = TRUE)+
  geom_text(aes(label=value_text) , position=position_stack(vjust =0.5) , size=4 , color="white")+
  xlab("") +
  scale_fill_manual(values=col) +
  theme(
    title =element_text(size=5, face='bold'),
    strip.text.x = element_text(size = 10, colour = "black",face='bold') ,
    axis.text.x = element_text(size = 8,color="black",face='bold') ,
    axis.text.y = element_text(size = 8,color="black",face='bold') ,
    axis.title.y = element_text(size = 10,color="black",face='bold') ,
    #axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1 , size = 8,color="black",face='bold'),
    legend.title = element_blank(),
    legend.text = element_text(size = 8),
    legend.position = "bottom",
    axis.title =element_text(size = 8),axis.text =element_text(size = 7, color = 'black')
  )


images_name <- paste0(out_path , "/Ratio_molecular_subtype.pdf")
ggsave( images_name , p , width = 8 , height = 5 )
